C -------------------------------------------------------------------------
C SUBROUTINE lorenz(x,y,dydx) RETURNS DERIVATIVES dydx AT x FOR THE
C LORENZ ATTRACTOR
C -------------------------------------------------------------------------
      
      subroutine lorenz (x,y,dydx)
C     =================

      implicit real*4 (a-h,o-z)
      parameter (nmax=10) 
      dimension y(nmax),dydx(nmax)

      common /parm/ sig,b,r

      dydx(1)=sig*(-y(1)+y(2))
      dydx(2)=-y(1)*y(3)+r*y(1)-y(2)
      dydx(3)=y(1)*y(2)-b*y(3)

      return
      end


C -------------------------------------------------------------------------
C GIVEN VALUES FOR n VARIABLES y AND THEIR DERIVATIVES dydx KNOWN AT x, USE 
C THE FOURTH ORDER RUNGE-KUTTA METHOD TO ADVANCE TO THE SOLUTION OVER AN
C INTERVAL h AND RETURN THE INCREMENTED VALUES AS yout, WHICH NEED NOT BE
C A DISTINCT ARRAY FROM y.
C SUBROUTINE derivs(x,y,dydx) RETURNS DERIVATIVES dydx AT x FOR n
C DIFFERENTIAL EQUATIONS
C -------------------------------------------------------------------------
      
      subroutine rk4 (y,dydx,n,x,h,yout)
C     ==============  

      implicit real*4 (a-h,o-z)
      parameter (nmax=10) 
      dimension y(n),dydx(n),yout(n),yt(nmax),dyt(nmax),dym(nmax) 

      hh=h*0.5
      h6=h/6. 
      xh=x+hh 

      do 11 i=1,n 
        yt(i)=y(i)+hh*dydx(i) 
11    continue

      call lorenz(xh,yt,dyt)

      do 12 i=1,n 
        yt(i)=y(i)+hh*dyt(i)
12    continue

      call lorenz(xh,yt,dym)

      do 13 i=1,n 
        yt(i)=y(i)+h*dym(i) 
        dym(i)=dyt(i)+dym(i)
13    continue

      call lorenz(x+h,yt,dyt) 

      do 14 i=1,n 
        yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
14    continue

      return
      end 


C -------------------------------------------------------------------------
C RUNGE-KUTTA QUALITY CONTROL
C FIFTH-ORDER RUNGE-KUTTA STEP WITH MONITORING OF LOCAL TRUNCATION ERROR
C TO ENSURE ACCURACY AND ADJUST STEPSIZE. INPUT ARE THE DEPENDENT VARIABLE
C VECTOR y OF LENGTH n AND ITS DERIVATIVE dydx AT THE STARTING VALUE OF THE
C INDEPENDENT VARIABLE x. ALSO INPUT ARE THE STEPSIZE TO BE ATTEMPTED htry
C THE REQUIRED ACCURACY eps, AND THE VECTOR yscal AGAINST WHICH THE ERROR
C IS SCALED. ON OUTPUT, y AND x ARE REPLACED BY THEIR NEW VALUES, hdid IS 
C THE STEPSIZE WHICH WAS ACTUALLY ACCOMPLISHED, AND hnext IS THE ESTIMATED
C NEXT STEPSIZE. 
C -------------------------------------------------------------------------

      subroutine rkqc (y,dydx,n,x,htry,eps,yscal,hdid,hnext)
C     ===============

      implicit real*4 (a-h,o-z)
      parameter (nmax=10,fcor=.0666666667,one=1.,safety=0.9,
     &           errcon=6.e-4) 

C errcon EQUALS  (4/safety)**(1/pgrow)

      dimension y(n),dydx(n),yscal(n),ytemp(nmax),ysav(nmax),dysav(nmax)

      pgrow=-0.20 
      pshrnk=-0.25
      xsav=x

      do 11 i=1,n 
        ysav(i)=y(i)
        dysav(i)=dydx(i)
11    continue

      h=htry
1     hh=0.5*h

      call rk4(ysav,dysav,n,xsav,hh,ytemp) 

      x=xsav+hh 
      call lorenz(x,ytemp,dydx) 
      call rk4(ytemp,dydx,n,x,hh,y)

      x=xsav+h
C      if(x.eq.xsav)pause 'stepsize not significant in rkqc.'

      call rk4(ysav,dysav,n,xsav,h,ytemp)

      errmax=0. 
      do 12 i=1,n 
        ytemp(i)=y(i)-ytemp(i)
        errmax=max(errmax,abs(ytemp(i)/yscal(i))) 
12    continue
      errmax=errmax/eps 

      if(errmax.gt.one) then
        h=safety*h*(errmax**pshrnk) 
        goto 1
      else
        hdid=h
        if(errmax.gt.errcon)then
          hnext=safety*h*(errmax**pgrow)
        else
          hnext=4.*h
        endif 
      endif 

      do 13 i=1,n 
        y(i)=y(i)+ytemp(i)*fcor 
13    continue

      return
      end 


C -------------------------------------------------------------------------
C RUNGE-KUTTA DRIVER WITH ADAPTIVE STEPSIZE CONTROL.
C INTEGRATE THE nvar STARTING VALUES ystart FROM xstart WITH ACCURACY
C eps TAKING nstep STEPS. 
C hdes IS THE DESIRED STEPSIZE FOR WHICH THE RESULTS ARE WRITTEN TO FILE,
C hmin IS THE MINIMUM ALLOWED STEPSIZE (CAN BE ZERO). 
C -------------------------------------------------------------------------

      subroutine odeint 
     & (ystart,nvar,xstart,nskip,nstep,path,eps,hdes,hmin, psig, pb, pr) 
C     =================

      implicit real*4 (a-h,o-z)

      parameter (nmax=10,maxstp=10000)
      parameter (two=2.0,zero=0.0,tiny=1.e-30) 

      dimension ystart(nvar),yscal(nmax),y(nmax),dydx(nmax)
      dimension path(nmax,maxstp)

      common /parm/ sig,b,r

      x=xstart
      kount=0

	sig = psig
	b = pb
	r = pr

      do 11 i=1,nvar
        y(i)=ystart(i)
11    continue

C LOOP OVER STEPS
      do 16 nstp=1,nskip+nstep

        h=hdes
        hsum=0.
20      continue

C SET SCALE USED TO MONITOR ACCURACY
        call lorenz(x,y,dydx) 
        do 13 i=1,nvar
          yscal(i)=abs(y(i))+abs(h*dydx(i))+tiny
13      continue

C CHECK COUNTER
        kount=kount+1
C        if(kount.gt.maxstp) pause 'max. number of steps exceeded.'

C GET NEW VALUES
        call rkqc(y,dydx,nvar,x,h,eps,yscal,hdid,hnext)
C        if(abs(hnext).lt.hmin) pause 'stepsize smaller than minimum.' 
        hsum=hsum+hdid

C ONE MORE STEP ?
        if(hsum.lt.hdes)then
          h=min((hdes-hsum),hnext)
          goto 20
        endif      

C STORE INTERMEDIATE RESULT
        if ( nstp.gt.nskip ) then
          do 15 i=1,nvar
            path(i,nstp-nskip)=y(i)
15        continue
        endif
  
16    continue

      return
      end 

